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We theoretically study the dynamic screening properties of bilayer graphene within the random phase ap- 
proximation assuming quadratic band dispersion and zero gap for the single-particle spectrum. We calculate the 
frequency dependent dielectric function of the system and obtain the low energy plasmon dispersion and broad- 
ening of the plasmon modes from the dielectric function. We also calculate the optical spectral weight (i.e. the 
dynamical structure factor) for the system. We find that the leading order long wavelength limit of the plasmon 
dispersion matches with the classical result for 2D electron gas. However, contrary to electron gas systems, the 
non-local plasmon dispersion corrections decrease the plasmon frequency. The non-local corrections are also 
different from the single layer graphene. Finally, we also compare our results with the double layer graphene 
system (i.e. a system of two independent graphene monolayers). 
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I. INTRODUCTION 

Bilayer graphene (BLG) has recently attracted experimen- 
tal and theoretical interest, both for fundamental physics as 
well as for possible technological applications' The system 
consists of two layers, each of which has carbon atoms ar- 
ranged in a honeycomb lattice with two (A and B) sublattices. 
Each layer can be viewed as a single layer graphene (SLG)^ 
with their associated physics of chiral Dirac Fermions, cou- 
pled by strong interlayer tunneling between A and B lattice 
sites in the two layers. It is useful to view BLG as an effective 
layer where the interlayer tunneling changes the linear disper- 
sion of the quasiparticles in SLG to an effective quadratic dis- 
persion at small wavevectors^. Thus, a simple model for BLG 
single-particle spectrum is a pair of chiral parabolic electron 
and hole bands touching each other at the Dirac (or the charge 
neutrality) point. Each band has a 4-fold degeneracy arising 
from spin and valley degrees of freedom. The pristine or un- 
doped bilayer graphene has attracted a lot of interest since the 
presence of a single Fermi point and quadratic dispersion can 
lead to a host of exotic phenomena^ ^. The low energy prop- 
erties of the doped bilayer graphene also show interesting fea- 
tures like enhanced backscattering due to the chirality of the 
bands present in the system"^. The single-particle spectrum 
is sufficiently different in SLG (linear dispersions) and BLG 
(quadratic dispersions) that it is theoretically interesting to ex- 
plore and contrast various dielectric properties in the two sys- 
tems. In this paper we consider the dynamic screening and 
low energy collective modes of BLG. 

Coulomb interaction and its dynamic screening due to 
many-body effects is an important property of any electronic 
material. While static screening determines transport proper- 
ties, the dynamic (frequency dependent) screening determines 
the elementary quasiparticle spectra and the collective modes 
and is crucial to understanding the optical properties of the 
system. While much work has been done on screeningprop- 
erties and plasmon modes of SLG both theoreticall}M] and 
experimentally'^, there are very few existing analytic works 
on doped bilayer systems. There are numerical calculations 
of the bilayer dynamic screenin g with in a four-band model'', 
and the parabolic approximation'- but it is not useful in a 



general context. Re f.T? obtains results which are qualitatively 
similar to the results in this paper. Recently, the static screen- 
ing properties of BLG systems have been studied analytically 
in Ref. 9 In this paper, we will analytically study the dynamic 
screening of Coulomb interactions in doped BLG and obtain 
its collective modes. 

The parabolic dispersion makes the BLG low energy 
physics quite different from SLG systems. On the other hand, 
although the system has the same low energy dispersion as a 
regular two dimensional electron gas (2DEG) system, the chi- 
rality of graphene leads to features which are distinct from a 
standard 2DEG'^ which has been extensively studied in the 
context of semiconductor heterostructureJ^. It is therefore 
useful to compare the physics of BLG with SLG and 2DEG. 
In this respect, it is important to note that although bilayer 
graphene consists of two layers of graphene sheets, it should 
really be viewed as a 2D system with parabolic dispersion and 
chiral bands for the purpose of low energy physics. Thus, it is 
also important to compare this sys tem w ith other bilayer sys- 
tems such as double quantum 

wellPSa 

which consists of two 
layers of 2DEG coupled by Coulomb interaction and shows 
distinct features as a result of its quasi-2D nature. Another 
closely related system is the double layer graphene, which is 
obtained by putting a layer of oxide between two single lay- 
ers of graphene'^. In this case, the layer separation is much 
larger (~ lOnm) than in BLG (^ O.Snm) and the system is 
better viewed as two different layers of graphene coupled by 
Coulomb interactions without any interlayer tunneling. The 
dynamic screening properties of BLG is quite different from 
these systems. It is then useful to compare the results obtained 
for BLG with corresponding results in these two bilayer sys- 
tems (i.e. two layers of 2DEG and two layers of SLG). 

In this paper, we analytically study the dynamic screen- 
ing properties of the Coulomb interaction in BLG systems 
within the random phase approximation (RPA). We calculate 
the dielectric function of bilayer graphene, e(q, w) at arbi- 
trary wavevectors q and frequency lo. The imaginary part 
of l/e(q, is related to the optical spectral weight which 
can be directly measured in optical conductivity or light- 
scattering^^, or electron-scattering measurements"^' while 
the zeroes of the real part give the dispersion of the plasmon 
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FIG. 1: (a)-(c): Imaginary part of Intraband susceptibility as a function of frequency for (a) q = O.S/cf (b) q — l.Skp and (c) q = 2.5fc_F. 
(d)-(f): Real part of Intraband susceptibility as a function of frequency for (d) q — O.Skp (e) g = LSfc^ and (f) q — 2.5A;_f. 



modes. We calculate the plasmon modes and their damping 
due to the presence of the second band in the system. We 
make extensive comparisons of our results with similar results 
for SLG systems'*^ and 2DEG'^, especially regarding the de- 
pendence of collective mode dispersions on the carrier densi- 
ties in these systems. We look at the density dependence of the 
long wavelength collective mode dispersion and show how the 
system crosses over from a BLG like behaviour to a SLG Uke 
behaviour as density of the carriers is increased. We also com- 
pare our results with those for quasi-2D systems like double 
quantum well^ and double layer graphene^^ and show that 
in spite of having two layers, bilayer graphene shows unique 
features in its dynamic screening properties. 

We note that we are working here with the parabolic ap- 
proximation to the dispersion of the BLG bands. In real bi- 
layer systems, the dispersion changes from parabolic to a lin- 
ear form (i.e. the full dispersion is hyperbolic) at high en- 
ergies. In addition the 2 band approximation breaks down 
at higher energies, resulting in a 4-band BLG model. How- 
ever, it is clear that for low density materials, the low energy 
properties of the system are qualitatively well described by 
keeping only the parabolic part of the dispersion. For bilayer 
graphene systems, our parabolic approximation for collective 
modes should be valid upto densities ^ 5 x lO^^cm^^. 

Our BLG model, as mentioned above, consists of a 2-band 
(single valence and conduction band) gapless parabolic chi- 
ral 2D single-particle energy dispersion. We compare our an- 
alytic dynamic screening and collective mode BLG results 
with SLG (chiral gapless 2-band linear dispersion), 2DEG 
(non-chiral 1 -band gapped parabolic dispersion), double-layer 
graphene (two parallel SLG) and 2D bilayer ( two parallel 
2DEG) systems. Our theoretical results, therefore, give a 



fairly complete picture of dielectric screening and plasmon 
modes in 2D systems, covering both single and bilayer sys- 
tems, both linear and quadratic band dispersions, system with 
and without gaps, and both chiral and non-chiral situations. 
Our calculated BLG plasmon dispersion can be directly com- 
pared with experimental results when they become available. 



The low energy effective Hamiltonian for BLG is given by 



H 



1 



(fcj; iky)^ 



2m \ {kx + ikyY 







(1) 



where ra = j/{2vp) , 7 is the interlayer tunneling amplitude 
inherent in the BLG system and vp is. the graphene Fermi ve- 
locity. We set h — 1 throughout this paper. Typically, in the 
BLG systems m ~ 0.033mo, with tuq being the free electron 
mass. This Hamiltonian differs from the SLG case in having 
a quadratic as opposed to a linear dispersion. As a result, the 
system has a constant density of states rather than the linear 
density of states in SLG. As we will see, this leads to substan- 
tial difference in the BLG dielectric response. 



The Hamiltonian can be diagonalized to obtain two bands 
with dispersions — sk^/(2TO) and coiTesponding wave- 
functions ^1= ^ (p-'^^^ 



^,s), where s = ±1 denotes re- 



spectively the conduction and the valence band. Here Xk — 
t&n^^ (ky/kx). In this paper we use the parabolic Hamilto- 
nian to obtain the dynamic screening within random phase 
approximation (RPA). This involves a theoretical calculation 
of the dynamical dielectric function e(q, uj). 
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FIG. 2: (a)-(c): Imaginary part of Interband susceptibility as a function of frequency for (a) q = O.S/cf (b) q — l.Skp and (c) q = 2.5fc_F. 
(d)-(f): Real part of Interband susceptibility as a function of frequency for (d) q — O.Skp (e) g = LSfc^ and (f) q — 2.5A;_f. 



II. BARE BUBBLE POLARIZABILITY 



The frequency dependent dielectric function for this system 
(within RPA) is given by 



e(q,w =1 n q,w 



(2) 



where k is the background dielectric constant. Here 11 is the 
free particle (bare bubble) polarizability, given by 

n(q,w)=.g^^^^^%-F,,,(k,k + q), (3) 

where is the Fermi function, g is the degeneracy factor 
(g — 4) and the wavefunction overlap Fgs' (k, k + q) = 
i[l + ss' cos{26)] with 6 being the angle between k and 
k + q. Using a circular co-ordinate system {k,(j>), one can 
write cos 9 = {k + (7cos0)/|k + q|. The chirality effect is 
incorporated in the matrix element F. 

The polarizability can be separated into the intmband is = 
s' terms) and the interband (s' = —s terms) contributions. 
The scale for the polarizability is set by the free particle den- 
sity of states A^o — gm/(27r). Using the Fermi momentum 
kp^ and the Fermi energy Ep as, units of length and en- 
ergy respectively, we work with the dimensionless variables 
X = k/kp, y = q/kp and z = ujjEp. Then, we can write 
n(g, w) = iVo[n™*'=''(y, z) + n™*™(?/, z)], where n™*™and 
jjmter dimcnsionlcss quantities. Note that the interband 
contribution is important in graphene only because of its gap- 
less nature; in semiconductor based 2DEG, the interband con- 



tribution is implicitly included in the background dielectric 
constant and is not considered in the electronic polarizability. 



We consider a gated or doped BLG with Ep ^ 0, where 
the Fermi level lies in the conduction band. At T = 0, this 
implies that = Q{kp ~ k) and — ^- Then one can 
write n"'*™(t/, z) = U^{y, z) + ni(-y, -z), where 



n^(?/, z) — — j xdx [ c 



1 



1 - 



z — y"^ ~ 2xy cos ( 



y"^ sin^ ( 



x'^ + y"^ + 2xy cos cf) 
and n"*'='^(y, z) = n'^{y, z) + W^i-y, -z), where 
-1 



(4) 



sm 



xdx / d(j)—^ , „ , 
^ + 2/^ + Zxy cos <p 



y 



z + 2x^ + + 2xy cos < 



(5) 



The azimuthal integrals in eqns. Q and (jSjl can be done ana- 
lytically to obtain 



2 1 I 2 2| 

X + z — \x — y \ 



(6) 



Sgn[z — + 2xy] 



(2a;2 + z - 2/2)2 
\/{z — y^Y' ~ 4x2y2 



4 



dx 



x{x'^ + z) 
Sgn[{x — yY + x^ + z] 
a/(2x2 + Z + y2)2 _ 4a;2y2 



(7) 



where Sgn{x) — ±1, depending on a; ^ 0. The real and 
imaginary parts of the polarizabihty are even and odd under 
z — > — z; so we will only consider the case of z > 0. The real 
and imaginary parts of the intraband polarizabihty are 



Re n™*™(y,z) 



Im W 



fi{y,z) + fiiy, -z) 

e[{z~yY-^y']f2{y,z) 

e[{z + yY~^yV2{y,~z), (8) 



\y,z) = -e[4y'-{z-y'f]My,z) 

+ e[4y^-{z + y^f]f3{y,-z). (9) 



For notational brevity, it is useful to define the quantities a 

{z — y'^)/2y and h = {z ^ y^) I2y. Then we have 



ya 
2z 
ya 



log 



1 - z 
e[l-y] log 



1 

' 2 
1 - z 



2ya 



(10) 



f2{y,z) = Sgn[a 



— 1 yb 



2z 



log 



Va^ -l + b 



ya 
2z 



log 



1 -a 



1 + a 



(11) 



My,z) 



\/l - a^ y\a\ 



tan 



z 



VT 



\b\ 



(12) 



For the interband polarizabihty, it is useful to define: c = 

y^/y^~+~2z and d ~ y\/2z — y"^. Then we can write 

Re n™*-(y, z) = U{y, z) + /4(y, -z) + ^(y, z) 

+ Q[{z-2f + d'']Sgn{2-z^\d\)h{y, 



where 



1 


log 


c 






2 




4 



Z+J^ 

2z 



log|l + z| 



Vc^ + z2 

2^ 



log 



+ z2 



e[y-i]^^log 

2z 



Vc^ + z^ + z 
1 + z 



(15) 



z + y2 



1 
2 

log 



log 



2 + Z+ v/(2 + z)2 + c2 



2z 



Vc2 +Z2(2 + Z) - Zv/(2 + z)2+c2 



Vc2 +Z2(2 + Z) + Zv/(2 + z)2+c2 



(16) 



2z 



— tan 



zv/d2 - (2 - z)2 



(17) 



Im 



"(?/,z) = e[2z-?;2]e[z-2 + d; 
e[z-2 

+e[2 + d-z]/6(y,z) 



(14) 



Eqn.s (8p7i are the main analytic results of this paper. These 
equations completely define the BLG RPA dynamic response 
(and consequently the collective plasmon modes) analytically 
within the 2-band low energy quadratic band dispersion ap- 
proximation. We first point out the essential features of the 
polarizabihty which differentiate it from the 2DEG and SLG 
systems. 

The intraband polarizabihty has an imaginary part for < 
z < y"^ + 2y for y < 2 and for y-^ — 2y < z < y-^ + 2y 
for y > 2, which defines the band of particle-hole contin- 
uum. For y < 1, there is a discontinuity in the derivative at 
z = 2y — y^, which is qualitatively similar to the 2DEG re- 
sults. For 2 > y > 1, an additional derivative discontinuity 
appears at z — y^, which is a result of the chiral nature of 
the wavefunctions. For y > 2, the feature at z = y^ persists, 
while z = 2y — ?/2 shifts to the lower edge of the particle-hole 
continuum. These features are shown in Fig.[TJa), (b) and (c) 
respectively. These features result in sharp features in the real 
part of the polarizabihty. However, unlike SLG systems, the 
polarizabihty does not diverge at the upper edge of the con- 
(ISiinuum. The singularity in this case is softened to a derivative 
^^j") discontinuity. The real part of the polarizabihty for different 
wavevectors is plotted in Fig.[T](d), (e) and (f) respectively. 

The imaginary part of the interband polarizabihty is non- 
zero for z > j/2 _ 2y + 2 for y < 2 and for z > y"^ /2fox y > 2 
which defines the interband continuum. For z > y'^ + 2y + 2, 
it has the simple form ^ y^/z. For y > 1, there is a derivative 
discontinuity at z = j/2 xhe real part of the polarizabihty 
shows a peak at the lower edge of the interparticle continuum. 
The imaginary and real part of the interband polarizabihty are 
plotted in Fig.jSJa), (b), (c) and Fig.|2](d), (e) and (f) respec- 
tively. 
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FIG. 3: Optical plasmon dispersion in bilayer grapliene systems for (a) — 14 (b) = 5.6 and (c) Vs = 4.3. The tliick red line is the 
plasmon dispersion, while the dashed lines show the upper edge of the intraband continuum and the lower edge of the interband continuum, 
(color online). 




III. PLASMONS AND OPTICAL SPECTRAL WEIGHT 

The dielectric function can be used to calculate the disper- 
sion and damping of plasmons, which are the long-wavelength 
longitudinal collective modes of the density fluctuations. The 
dispersion can be obtained by looking for poles of the two- 
particle Green functions or equivalently for zeroes of the real 
part of the dynamical dielectric function. It is useful to express 
the dielectric function in terms of the dimensionless quantities 
y = q/kp and z =^ uj/Ep. 

e{q,Lu) - 1- ^(n»*™(y,z) +ff"*^'-(y,z)), (18) 

y 

where the dimensionless ratio = gm/ [nkp) denotes the 
ratio of the interaction energy scale to the kinetic energy scale 
in the system. It is important to note that unlike SLG systems, 
where is only a function of the substrate dielectric constant 
and can take a few discrete values, for BLG systems can be 
varied continuously by changing the gate voltages and hence 
density of the carriers. The coupling constant in BLG is 
thus similar to that in 2DEG with carrier density being the 
tuning parameter (r^ ~ n^^^^). 

There are two collective modes of this system, an optical 
plasmon mode which lies above the intraband electron-hole 
continuum and an acoustic plasmon mode which lies within 
the intraband continuum. The acoustic mode, which has a lin- 



ear dispersion at long wavelength limit is always overdamped. 
Thus it caiTies very little spectral weight and is not experimen- 
tally relevant. In the following we will focus on the optical 
plasmon mode (henceforth referred to as the plasmon). 

To get the dispersion of the optical plasmon, we note that 
the mode lies between the intraband and interband continuum 
in the low wavelength limit. In this limit we can expand the 
polarizability to obtain 

2 2 

n(y,z) = 2^-^. (19) 

In the long wavelength limit, we find that the plasmon has a 
leading order ^ y/q dispersion. In this limit, we recover the 
2DEG result in leading order in q 

Expanding to next order, we can obtain the quantum (non- 
local) corrections to the plasmon dispersion. The BLG plas- 
mon dispersion is upto this order is given by 

IgEpq ( qTFq\ 

where qTP = rgkp is independent of the density. 

We next compare the plasmon dispersion for BLG systems 
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FIG. 5: Plasmon frequency for BLG (solid black) and SLG (dashed 
red) systems for two different values of substrate dielectric constant 
n — 1 and K — 10. The two figures are for different densities (a) 
n = lO^^cm^^ and (b) n = 10^'^ cm^^ (color online). 



with those of SLG and a 2DEG. For the sake of comparison, 
we first note that the plasmon dispersion formulae for both 
SLG"* and 2DEG'^ can be cast in a form similar to BLG, 
cjq — ujoq^^^{l — aq/k^). We first focus on the leading order 
^ y/q term. For BLG , the coefficient of the leading order 
term is given by ujq — {ge^Ep/ kY^"^, while for SLG it is 
given by = {ge^vphp/'^Ky^^ . For the 2DEG it is given 
by ujq — (jje^Ep I k)^!'^ . The density dependence (~ n^/^) 



of Wo in BLG is different from that of (~ n^l\ and 

is a consequence of a constant as opposed to a linear density 
of states. The long wavelength BLG plasmon is thus identi- 
cal to the ordinary 2D plasmon and is thus different from the 
long wavelength SLG plasmon. In particular, the long wave- 
length BLG plasma frequency, in contrast to the SLG plasma 
frequencj'^, is classical and does not have any h in the leading 
order dispersion. 

However, the difference with the 2DEG result occurs in 
the next order which encapsulates the effect of the quan- 
tum corrections. For the BLG plasmon, a = qTp/8 with 
qpp = rgkp being independent of density. For the 2DEG, 
a = SqTp/-i with qxp = grne^ / k. Thus, the quantum 
corrections increase the plasmon frequency from its leading 
order result in 2DEG, while in BLG they lead to a reduc- 
tion of the plasmon frequency. For SLG, a = qxp/^ with 



Qtf = ge^kp/vpK ~ n^/^. The quantum corrections thus 
decrease the SLG plasmon frequency. However the density 
dependence of this correction term is very different from the 
BLG case. 

The other major difference from the single layer graphene 
is that in SLG systems, there is a singularity of the real part 
of the polarizability at the upper edge of the intraband con- 
tinuum. Hence, the plasmon dispersion never goes into the 
continuum and the pole survives upto arbitrary large wave- 
vectors. In bilayer systems, there is no discontinuity at the 
upper edge of the intraband continuum (it is replaced by a dis- 
continuity in the derivative with respect to frequency), and so, 
at a critical wave-vector, the plasmon pole vanishes once the 
plasmon dispersion touches the edge of the continuum. 

In order to compare BLG and SLG plasmon dispersions 
we plot in Fig.[5]the plasmon dispersion (in absolute units of 
eV) as a function of wavevector (in absolute units of nm^^) 
for two fixed density of careers (a) n = lO^'^cm^^ and (b) 
n = 10^"^ cm" ^ . For each value of density we plot the plas- 
mon dispersion for two values of k = 1 and k — 10. We see 
that the BLG plasmon merges with the intraband continuum 
at a finite wavevector and gets overdamped while the SLG 
plasmon is present at all wavevectors. At the lower density 
the SLG plasmon frequency is higher than the BLG plasmon 
while at the higher density the BLG plasmon has a higher fre- 
quency than SLG plasmon as wq oc ^/Ep and the Fermi en- 
ergy has different dependence on density in the two cases. 

At long wavelengths, the plasmon mode is completely un- 
damped. As the dispersion goes into the interband continuum 
at larger wavevectors, two things happen: (i) the mode gets 
damped due to usual Landau damping by the particle-hole 
pairs, and (ii) the dispersion develops a knee like structure 
due to mode repulsion from the continuum. 

In order to look at typical BLG plasmon dispersions, we 
consider three types of substrates: (a) SiC, which has k = 
5.5, (b) Si02, which has k — 4 and (c) vacuum, correspond- 
ing to suspended bilayer graphene, with k = 1. Note that 
the effective k to be used for calculation of is the average 
of K for the substrate and that for the vacuum. For a typical 
density of lO^^cm"^, we get ^ 4.3 for SiC, = 5.6 for 
Si02 and r., = 14 for suspended bilayer graphene. The plas- 
mon dispersion for these three values of are shown in Fig. [3] 
The critical wavevector for plasmon damping does not change 
much varying from qi = 0.29kp for — 14 to qi = 0.38kp 
for Ts — 4.3. The maximum wavevector for which the plas- 
mon pole exists varies from qmax = l-12kp for = 4.3 to 
qmax = l-52kp for rs = 14. 

An important quantity which is directly measured in spec- 
troscopic experiments is the optical spectral weight or the loss 
function. This quantity is proportional to the imaginary part 
of the inverse dielectric function of the system. In Fig. [4j we 
plot the spectral weight for three different wavevectors for 
rs = 4.3. The first wavevector {q ~ Q.15kp) coiTesponds 
to the situation where the plasmon mode lies between the in- 
traband and interband continuum. The plasmon in this case 
is sharp, i.e. not Landau-damped. The second wavevector 
(q — O.Skp) corresponds to the situation where the plas- 
mon peak has entered into the interband continuum. In this 
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FIG. 6: (a) The quasiparticle residue Z and (b) The width F of the 
plasmon pole in the inverse dielectric function for ~ 4.3. 
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FIG. 7: The scaling of the leading order plasmon dispersion in the 
long wavelength limit with density on a log-log scale. The leading 
order result is obtained by evaluating the plasmon dispersion at a 
wavevector 5 = 5 x lO^cm^^. The red line is the BLG result with 
parabolic dispersion. The blue line is the SLG result. The black 
line is the BLG plasmon dispersion with hyperbolic dispersion (color 
online). 



can be mapped out as a function of q and oj. 



case we see a broadened plasmon peak. The third wavevector 
(q = 2.2kp) corresponds to the case where the plasmon pole 
has vanished, and the collective mode is overdamped. It is to 
be noted that since the dielectric function continues to have 
a kink at the edge of the intraband continuum even when the 
plasmon has ceased to exist, a peak like structure is seen in 
the optical spectral weight in this limit. This peak vanishes at 
much larger wavevectors ^ Akp for — 4.3. 

Focusing on the plasmons, one can expand the inverse di- 
electric function near the plasmon pole to write 



1 



UJ — UJg + iTq 



(22) 



where the 



is the plasmon frequency, Zq 



~l/[V{q)dIl {q,LL!)/duj\uj=ujJ is the plasmon pole residue or 
equivalently, the weight of the plasmon peak, and the width of 
the peak is given by Fq = — 1/(17)11 (q, ajg)Zq. Here H and 
n denote the real and imaginary part of the polarizability. In 
Fig. [6] we plot Zq and F^ as a function of q for = 4.3. 

The plasmon weight Zq ^ ^/q in the long wavelength limit. 
It shows a broad peak before the plasmon dispersion enters the 
interband continuum. At the lower edge of the interband con- 
tinuum, Zq shows a kink as a function of q. It flattens out 
over a substantial range of wavevectors before finally going 
to zero as the plasmon approaches the intraband continuum 
and vanishes. The width of the plasmon peak starts from zero 
at the interband edge and shows a broad feature before going 
to zero again as the plasmon vanishes. The plasmon features 
predicted in our theory (Figs. |3]-|6]l should be directly observ- 
able in inelastic electron energy loss (lEEL)"^' or inelastic 
light scattering (ILS) spectroscopies^^ where the loss function 



IV. HYPERBOLIC DISPERSION AND CROSSOVER 
FROM BLG TO SLG PHYSICS 

We have treated bilayer graphene as a system with quadratic 
band dispersion and chiral electron and hole bands. While the 
parabolic dispersion is relevant at low energies, the actual dis- 
persion of the bands is hyperbolic. Thus at large wavevectors 
(relevant at large densities) the dispersion is effectively lin- 
ear and the system should behave like single layer graphene, 
while at small wavevectors (relevant at low densities) it should 
behave like the bilayer graphene system we have described. In 
this section, we wish to study this crossover by focusing on the 
long-wavelength limit of the plasmon dispersion. Although 
the full plasmon dispersion at all wavevectors can be obtained 
numerically, working with the leading order long wavelength 
dispersion lets us obtain analytic results and hence provides 
insight into the crossover. This will also let us estimate the 
regime of validity of our parabolic band approximation. 

To estimate the long wavelength plasmon, we note that for 
qvp/i^ ^ 1, the bare bubble polarizability is given by 



n(q,L.) 



gkp q^ dEi, 



47r w2 



(23) 



k=kf 



where is the dispersion of the band containing the Fermi 
level. Hence the long wavelength plasmon dispersion is given 
by 



ge^kpqdEF ge^ndEp 



2 k dkp 



dn 



(24) 



where n is the density of carriers and we have used n ~ 
gkp/{4:7r). The above results are true for any dispersion in- 
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eluding (a) a parabolic dispersion, (b) a linear dispersion and 
(c) a hyperbolic dispersion. Let us now focus on the hyper- 
bolic dispersion of BLG. In this case one can write the Fermi 
energy as a function of density aP 



T 
2 



1/2 



(25) 



where 7 is the interlayer tunneling and vp is, the Fermi ve- 
locity of the graphene sheets. Note that in the limit of small 
n 0, Ep ^ kp/2m with m — ^/{2vp) and we re- 
cover the parabolic Fermi energy, while in the large density 
limit (n — )■ 00), Ep — vpkp and we recover the single 
layer graphene results. Now the long wavelength plasmon fre- 
quency for hyperbolic dispersion is then given by 



2 27rne^ Vp 



Ef 



1 - 



57^ 



(26) 



In the low density limit (n — > 0) this reduces to oj^ = 
{2Tine'^)q/ [nm) and we recover the leading order result for 
the bilayer graphene with parabolic dispersion. In the large 
density limit (n — > 00), the limiting behaviour is given by 
cjq = e^vp^gimq/ K which matches with the result for sin- 
gle layer graphene^. In between the answer smoothly inter- 
polates between these two limits. 

In Fig. |7] we plot the plasmon frequency at a low wavevec- 



tor {q = bx 10'' 



) as a function of density on a log-log 



scale. The low wavevector ensures that we are looking at the 
leading order result. We plot the results for BLG and SLG 
which show the usual n^/^ and m}/^ scaling respectively. The 
curve for hyperbolic dispersion of BLG smoothly interpolates 
between the above two curves. From the figure we can esti- 
mate the value of critical career density where the hyperbolic 
result differs significantly from the parabolic result. By com- 
paring the long wavelenght plasma frequencies in SLG and 
BLG we obtain the crossover density as ric — rn^Vp/TT. This 
correspond to a density n ^ 3 x 10^^cm~^ for a bilayer ef- 
fective mass m ~ 0.033me (see Fig. |7]l. This then gives us 
the limit of validity of the parabolic approximation for BLG 
systems. 



V. COMPARISON WITH OTHER BILAYER SYSTEMS 

The study of dielectric functions and plasmons for bilayer 
systems have a long history starting with the prediction of an 
undamped acoustic plasmon mode in double quantum well^ 
with layer separation greater than a critical distance. Simi- 
lar results have also been obtained for double layer graphene 
systems^"*. However, the physics of these systems is quite dif- 
ferent from the bilayer graphene system which is of interest 
in this paper. These other systems are comprised of well sep- 
arated 2D layers coupled only by the Coulomb interaction. 
Interlayer tunneling is absent in these systems. On the other 
hand, the bilayer graphene system is in the opposite limit of 
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FIG. 8: Low energy plasmon dispersion for three different bilayer 
systems: (a) an n-GaAs double quantum well with Ts = 1 (corre- 
sponding to a density of n = 3.2 x lO^^cm^^ and a layer separation 
of d = 9.5nm, which corresponds to the effective Bohr radius of 
n-GaAs quantum well) (b) a double layer graphene with a density 
of n = lO^^cm"^ and a layer separation of d = 3A (which corre- 
sponds to the layer separation of bilayer graphene) and (c) a bilayer 
graphene system with a density of n = 10^^cm~^ and = 4.3. 
The sqrtq mode is in blue and the linear mode is in red solid lines. 
The black lines indicate interband and intraband continuum. The un- 
damped linear acoustic mode is absent in bilayer graphene ( color 
online). 



small layer separation with layers coupled strongly through 
interlayer tunneling, which is accounted for in the model by 
changing the low energy dispersion from linear to quadratic. 
This leads to all the major differences between the bilayer 
graphene and other double layer systems. 

The double layer 2DEG systems in the limit of infinite sep- 
aration corresponds to two copies of the underlying 2D sys- 
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terns. Hence it has two plasmon modes (one in each layer), 
both of which have ^ ^ dispersion. For long wavelengths, 
these modes are undamped as they lie above the particle-hole 
continuum. The finite separation and the resultant Coulomb 
repulsion change the dispersion of one of the modes to lin- 
ear in q, but with a slope which keeps the plasmon above the 
two particle continuum. This is the origin of the undamped 
plasmon mode with linear dispersion'^. The slope of the lin- 
ear mode decreases with distance, and below a critical layer 
separation, the mode gets into the continuum and becomes 
damped. In double layer graphene we have the additional fea- 
ture that the undamped linear mode overlaps with the optical 
mode for a range of q values before splitting off as it enters 
the interband continuunH 

The crucial missing ingredient in the picture above is the 
effect of very strong interlayer tunneling which dominates the 
physics of bilayer graphene systems. In the presence of strong 
interlayer tunneling, the linear acoustic mode of the interac- 
tion coupled bilayer 2D systems is gapped out in BLG. In 
fact, the two band representation of BLG systems is an ef- 
fective low energy description. A more detailed description'^ 
would include 4 bands, two of which are gapped out (with 
gaps ^ eV). In this case one would recover the gapped acous- 
tic plasmon mode sitting at an energy above the highest band 
continuum. It is to be noted that the energy of the gapped 
mode is ^ eV and therefore does not affect the low energy 
properties of the system. Our low energy description of BLG 
system then recovers the low energy optical plasmon mode 
with ^ ^q dispersion. 

In Fig. [8] we show in the same figure, for the sake of com- 
parison, the low energy plasmon dispersion for typical param- 
eter values for the three bilayer systems. The double 2DEG 
dispersion is calculated for large layer separations and clearly 
shows the undamped optical and acoustic plasmon modes. 
The double layer graphene calculation is done for a relatively 
shorter layer separation. In this case the acoustic plasmon 
mode is quite close to the particle-hole continuum. Finally, the 
linear acoustic mode is not seen in bilayer graphene systems 
since the BLG, being strongly tunnel-coupled is effectively a 
1 -component system similar to a single 2DEG or SLG, but 
with its own distinct plasmon dispersion. 



VI. CONCLUSIONS 

In this paper, we have analytically studied the dynamic 
screening properties of doped bilayer graphene systems 
within RPA. We have obtained analytic forms for the full 
wavevector and frequency dependent dielectric function for 
the system at arbitrary densities. From this we have obtained 
the plasmon dispersion and calculated the plasmon spectral 
weight, focusing on the plasmon residue and broadening. We 
find that while the leading order plasmon dispersion matches 
with the classical 2D electron gas result, the non-local dis- 
persion correction in the next order suppresses the plasmon 
frequency, contrary to the 2D electron gas. However, the den- 
sity dependence of the Thomas Fermi wavevector is different 
from a single layer graphene, which is a consequence of the 



constant density of states as opposed to linear density of states 
in single layer graphene. This leads to a n}/^ density depen- 
dence of the BLG plasma frequency in contrast to the n^^^ 
dependence for SLG. 

We have also compared the results for bilayer graphene sys- 
tems with those in double quantum wells and double layer 
graphene systems. We find that while the Coulomb coupling 
between the layers results in an undamped linearly dispersive 
acoustic plasmon mode in double quantum wells and double 
layer graphene systems, the interlayer tunneling, which domi- 
nates the physics in the bilayer graphene systems, gap out the 
acoustic mode and we are left with a single undamped optical 
plasmon mode (dispersing as y/q) in the low energy limit. 

In conclusion, we discuss the approximations and the lim- 
itations of our theory, and in the process point out possible 
future research directions if experimental data on BLG plas- 
mon spectra, when they become available, necessitate further 
extension of our theory. We have used a T = RPA many- 
body theory for obtaining BLG dynamic dielectric response 
and calculating BLG plasmon spectra (plasmon dispersion, 
damping, and spectral weight). We have also assumed a clean 
system, neglecting any disorder effect on the dielectric re- 
sponse and the plasmons. In addition, we have made a 2-band 
parabolic single-particle band dispersion approximation. All 
of these approximations can, in principle, be improved albeit 
at the considerable cost of losing the analytic simplicity of our 
current theory. 

Given that the BLG Fermi temperature for a given carrier 
density n is given by Tp ~ AOOnK, where ft ~ n/lQ^^cm^^, 
our T = theory should be quite good even at room temper- 
ature for n > 10^^cm~^ and down to n ~ 10^^cm~^ (or 
perhaps even lower) as long as T > iK. The disorder effect 
is fairly easily incorporated in our theory (at least in the lead- 
ing order) by interpreting the frequency uj to be modified by 
the substitution cj^ — > + iTi), where = h/ {2Ti) is ob- 
tained from the system mobility /i (= a/ (ne), where a is the 
conductivity) by writing /i = er^ /m. The main effect of disor- 
der is thus to introduce a plasmon damping of Ti even outside 
the electron-hole Landau continuum regime. In the presence 
of Landau damping of plasmons, the disorder-induced adds 
to the plasmon broadening arising from Landau damping. In 
absence of Landau damping and at low temperatures, the dis- 
order broadening F^ is the main mechanism contributing to 
plasmon damping. 

In contrast to thermal and disorder effects, our other ap- 
proximations (i.e. RPA and parabolic dispersion) are not easy 
to improve upon. The parabolic band structure approxima- 
tion of our theory would fail at high energy where the BLG 
single-particle dispersion becomes linear, similar to SLG. 
This can be incorporated numerically in the RPA dielectric 
function, thus sacrificing the analytic transparency of our re- 
sults. We refrain from doing so because the SLG plasmon 
dispersion and dielectric response have already been calcu- 
lated in the literature When the BLG Fermi level is high in 
the conduction or valence band, the BLG dynamical response 
and plasmon mode dispersion will smoothly crossover to the 
SLG result. This is expected to happen for a carrier density 
> 5 X lO'^^cTO"^, where Ep is large enough for the BLG 
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band dispersion to be better approximated by a linear than a 
parabolic dispersion. 

Finally, RPA ignores many-body corrections such as vertex 
and self-energy effects in the polarizability, and includes only 
the important long-range divergence of the Coulomb interac- 
tion through the infinite series of bubble (or ring) polarizabil- 
ity diagrams. While RPA becomes exact for small q/kp (for 
all Ts) and for all q/kp (for rg <C 1), higher order self-energy 
and vertex corrections should become important for the po- 
larizability function n(q, lS) for large q/kp and large Vg- Un- 
fortunately, no systematic improvement of RPA dynamical di- 
electric response is available for arbitrary Vs values since such 
a theory must somehow include consistent and conserving in- 
finite series of many body diagrams. In the literature, unsatis- 
factory approaches based on "local field" corrections (i.e. cor- 
rections associated with correlation effects at large wavevec- 
tors transcending the zero wavevector Coulomb divergence) 
are often used to go beyond RPA, but theories involving lo- 
cal field corrections are not consistent from a diagrammatic 
many-body viewpoint since many-body diagrams from differ- 
ent perturbative orders are mixed in an uncontrollable manner. 
We can easily incorporate the most popular form of the local 
field approximation, the so-called Hubbard approximatioJ^, 



by rewriting eqn.|2]as 

e(q,w) = 1 - 



n(q,^) 



{l + ^G(<z)n(q,.;)} 



Kq 



(27) 



where the term in the denominator in curly bracket is the 
local field correction due to many-body effects neglected in 
RPA. The explicit form for the actual local field term G{q) 
is somewhat arbitrary except that G{q) — for q — Q, 
so that the RPA is recovered in the q limit. In 2D, 
G{q) = (l/2)(g/A/g2 + + q?^p) is often used although 
the rigorous validity of such a local field term in improving 
RPA is unknown. We do not pursue this line of reasoning be- 
cause such local field corrections are uncontrolled and static 
so that it is completely unclear whether it is an improvement 
on the dynamical RPA. If future BLG plasmon experiments 
show systematic deviations from our RPA predictions with in- 
creasing Tg (i.e. decreasing density, where our parabolic band 
approximation becomes better), then going beyond RPA may 
become necessary. 
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